####################################################################################
# Elections Improve Support for State Trial Court Judges in the United States ######
# 2023 Survey Analyses #############################################################
# Andrew R. Stone and Michael P. Olson #############################################
####################################################################################

####################################################################################
# Set working directory, load packages, functions ##################################
####################################################################################
  # Set working directory
  # Load packages
  library(cjoint); library(ggplot2); library(interflex); library(questionr); library(sandwich); library(stargazer); library(tidyverse); library(vtable)

  # Prediction and robust SE function
  predict.rob <- function(x,clcov,newdata){
    if(missing(newdata)){ newdata <- x$model }
    tt <- terms(x)
    Terms <- delete.response(tt)
    m.mat <- model.matrix(Terms,data=newdata)
    m.coef <- x$coef
    fit <- as.vector(m.mat %*% x$coef)
    se.fit <- sqrt(diag(m.mat%*%clcov%*%t(m.mat)))
    return(list(fit=fit,se.fit=se.fit))}
  
  # Ceiling function for p values
  p_ceiling <- function(x){return(ceiling(100*x)/100)}

#############
# Load data #
#############
  load("weidenbaum_survey_2023.RData")

##########################################################################
# Descriptive Statistics #################################################
##########################################################################
###########
# Table 2 #
###########
# Summary statistics for mechanism tests
  # Internal efficacy (outputting mean and SD)
  cat(as.character(format(round(mean(survey2023$cant_understand_politics,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/internal_mean.tex")
  cat(as.character(format(round(sd(survey2023$cant_understand_politics,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/internal_sd.tex")
  # External efficacy
  cat(as.character(format(round(mean(survey2023$officials_dont_care,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/external_mean.tex")
  cat(as.character(format(round(sd(survey2023$officials_dont_care,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/external_sd.tex")  
  
  # Elections force politicians to care
  cat(as.character(format(round(mean(survey2023$elections_make_govt_care,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/force_mean.tex")
  cat(as.character(format(round(sd(survey2023$elections_make_govt_care,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/force_sd.tex")
  # Local elections worth bothering with
  cat(as.character(format(round(mean(survey2023$local_elections_not_worth,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/local_mean.tex")
  cat(as.character(format(round(sd(survey2023$local_elections_not_worth,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/local_sd.tex")
  # Democracy satisfaction
  cat(as.character(format(round(mean(survey2023$democracy_satisfaction,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/satisfaction_mean.tex")
  cat(as.character(format(round(sd(survey2023$democracy_satisfaction,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/satisfaction_sd.tex")
  
  # Experience with democracy
  # This calculation is done in the real_world_experience_analysis.R file
  
  # Ideological distance
  cat(as.character(format(round(mean(survey2023$ideology_distance_5,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/ideology_mean.tex")
  cat(as.character(format(round(sd(survey2023$ideology_distance_5,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/ideology_sd.tex")
  # Punitiveness distance
  cat(as.character(format(round(mean(survey2023$punitiveness_distance,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/punitiveness_mean.tex")
  cat(as.character(format(round(sd(survey2023$punitiveness_distance,na.rm=T),2),nsmall=2)),sep = '\n', file = "./figures_tables/punitiveness_sd.tex")
  
#############
# Table C.1 #
#############
# Summary statistics table for respondent characteristics

  survey2023$Experience <- recode(survey2023$elected_any_court,"1"="Any Court Elected in State","0"="No Courts Elected in State")
  
  survey2023$GENDER <- dplyr::recode(survey2023$GENDER,"1"="Male","2"="Female")
  survey2023$Sex <- factor(survey2023$GENDER,levels=c("Female","Male"))
  
  survey2023$RACETHNICITY <- recode(survey2023$RACETHNICITY,"1"="White","2"="Black",
                                "3"="Other/Multiple","4"="Hispanic","5"="Other/Multiple",
                                "6"="Asian")
  survey2023$`Race or Ethnicity` <- factor(survey2023$RACETHNICITY)
  
  survey2023$EDUC5 <- recode(survey2023$EDUC5,"1"="Less than High School","2"="High School or Equivalent",
                         "3"="Associates Degree or Some College","4"="Bachelor's Degree","5"="Post-Graduate Study or Professional Degree")
  survey2023$Education <- factor(survey2023$EDUC5,
                             levels=c("Less than High School",
                                      "High School or Equivalent",
                                      "Associates Degree or Some College",
                                      "Bachelor's Degree",
                                      "Post-Graduate Study or Professional Degree"))
  
  survey2023$INCOME4 <- recode(survey2023$INCOME4,"1"="Less than $30,000","2"="$30,000 to $60,000",
                           "3"="$60,000 to $100,000","4"="$100,000 or More")
  survey2023$Income <- factor(survey2023$INCOME4,
                          levels=c("Less than $30,000",
                                   "$30,000 to $60,000",
                                   "$60,000 to $100,000",
                                   "$100,000 or More"))
  
  survey2023$IDEO <- recode(survey2023$IDEO,"1"="Very Liberal","2"="Somewhat Liberal",
                        "3"="Moderate","4"="Somewhat Conservative","5"="Very Conservative")
  survey2023$Ideology <- factor(survey2023$IDEO,
                            levels=c("Very Liberal",
                                     "Somewhat Liberal",
                                     "Moderate",
                                     "Somewhat Conservative",
                                     "Very Conservative"))
  
  sumtab <- sumtable(data=survey2023,
                     vars=c("Sex","Race or Ethnicity","Party",
                            "Education","Income","Ideology","Experience"),
                     out="return",
                     title="Descriptive Statistics: July 2023 NORC Survey")
  
  sumtab <- stargazer(sumtab,summary=F,
                      title="Descriptive Statistics: July 2023 NORC Survey",
                      label="norc-demographics-2023",rownames = F)
  
  sumtab <- gsub("ccc","lcc",sumtab,fixed=T)
  sumtab <- gsub("Sex","\\\\textbf{Sex}",sumtab)
  sumtab <- gsub("Race or Ethnicity","\\\\textbf{Race or Ethnicity}",sumtab)
  sumtab <- gsub("Party","\\\\textbf{Party}",sumtab)
  sumtab <- gsub("Education","\\\\textbf{Education}",sumtab)
  sumtab <- gsub("Income","\\\\textbf{Income}",sumtab)
  sumtab <- gsub("Ideology","\\\\textbf{Ideology}",sumtab)
  sumtab <- gsub("Experience","\\\\textbf{Experience}",sumtab)
  # Saving the table
  cat(sumtab, sep = '\n', file = "./figures_tables/studytwosumtab.tex")
  
##########################################################################
# Treatment Balance ######################################################
##########################################################################  
#############
# Table D.1 #
#############  
  # Using OLS, regressing respondent characteristics on profile attributes
  # Limiting to respondents who got treatment assignments
  survey.name <- survey2023[!is.na(survey2023$judge_name),]
  
  black_outcome <- lm(black ~ judge_name + judge_race 
                      + election_treatment_binary + judge_margin 
                      + judge_tenure + judge_party 
                      + judge_sentencing, 
                      data=survey.name)
  
  hispanic_outcome <- lm(hispanic ~ judge_name + judge_race 
                         + election_treatment_binary + judge_margin 
                         + judge_tenure + judge_party 
                         + judge_sentencing, 
                         data=survey.name)
  
  college_outcome <- lm(college ~ judge_name + judge_race 
                        + election_treatment_binary + judge_margin 
                        + judge_tenure + judge_party 
                        + judge_sentencing, 
                        data=survey.name)
  
  married_outcome <- lm(married ~ judge_name + judge_race 
                        + election_treatment_binary + judge_margin 
                        + judge_tenure + judge_party 
                        + judge_sentencing, 
                        data=survey.name)
  
  employed_outcome <- lm(employed ~ judge_name + judge_race 
                         + election_treatment_binary + judge_margin 
                         + judge_tenure + judge_party 
                         + judge_sentencing, 
                         data=survey.name)
  
  income_outcome <- lm(INCOME ~ judge_name + judge_race 
                       + election_treatment_binary + judge_margin 
                       + judge_tenure + judge_party 
                       + judge_sentencing, 
                       data=survey.name)
  
  internet_outcome <- lm(INTERNET ~ judge_name + judge_race 
                         + election_treatment_binary + judge_margin 
                         + judge_tenure + judge_party 
                         + judge_sentencing, 
                         data=survey.name)
  
  democrat_outcome <- lm(democrat ~ judge_name + judge_race 
                         + election_treatment_binary + judge_margin 
                         + judge_tenure + judge_party 
                         + judge_sentencing, 
                         data=survey.name)
  
  republican_outcome <- lm(republican ~ judge_name + judge_race 
                           + election_treatment_binary + judge_margin 
                           + judge_tenure + judge_party 
                           + judge_sentencing, 
                           data=survey.name)
  # List of models
  balance_list <- list(black_outcome,hispanic_outcome,college_outcome,married_outcome,
                       employed_outcome,income_outcome,internet_outcome,democrat_outcome,
                       republican_outcome)
  # Stargazer table of output
  balance_sg <- stargazer(balance_list,
                          se=list(sqrt(diag(vcovHC(black_outcome,type="HC2"))), # HC2 matches the output of amce (the cjoint function)
                                  sqrt(diag(vcovHC(hispanic_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(college_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(married_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(employed_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(income_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(internet_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(democrat_outcome,type="HC2"))),
                                  sqrt(diag(vcovHC(republican_outcome,type="HC2")))
                          ),
                          dep.var.labels = c("Black","Hispanic","College","Married","Employed",
                                             "Income","Internet","Democrat","Republican"),
                          covariate.labels = c("Woman","White","Hispanic","Appointed",
                                               "60\\% Margin","80\\% Margin",
                                               "5 Year Tenure","15 Year Tenure",
                                               "Democrat","Not Shown","6 Year Avg. Sentence",
                                               "9 Year Avg. Sentence","Constant"),
                          keep.stat=c("n","adj.rsq"),
                          notes.append = FALSE,
                          star.char="*",star.cutoffs=0.05,notes.label = "",
                          font.size="scriptsize",column.sep.width="0cm",no.space=T,
                          notes="\\parbox[t]{0.95\\textwidth}{\\scriptsize \\textit{Note}: Estimates are from OLS regressions with robust standard errors
                        in parentheses. The unit of observation is the respondent. All outcomes are binary except 
                        income, which is measured on an 18-point scale and is treated as a continuous variable. $^{*}p<0.05$.}",
                          title="Balance Test Results: 2023 Survey",
                          label="balance2023"
  )
  # Saving the table
  cat(balance_sg, sep = '\n', file = "./figures_tables/balance2023.tex")
  
####################################################################################
# Conjoint analyses ################################################################
####################################################################################
#######################  
# Set up the conjoint #
#######################
  # Names of attributes to plot
  attribute.names.full <- c("Selection Method",
                            "Gender","Race/Ethnicity",
                            "Selection Margin","Tenure",
                            "Partisanship","Sentencing Behavior")  
  
  # Names of varying attributes to plot
  levels.test.full<-list()
  levels.test.full[["a.election_treatment"]]<-c("Appointed","Partisan election","Nonpartisan election") # Baseline category is set below
  levels.test.full[["b.judge_name"]]<-c("Man","Woman")
  levels.test.full[["c.judge_race"]]<-c("Black","White","Hispanic")
  levels.test.full[["d.judge_margin"]]<-c("52","60","80")
  levels.test.full[["e.judge_tenure"]]<-c("1","5","15")
  levels.test.full[["f.judge_party"]]<-c("Republican","Democrat","Not Listed")
  levels.test.full[["g.judge_sentencing"]]<-c("Average","Below","Above") # Baseline category is set below
  
  # Adjusting baselines so sentencing baseline is 6 years (the average) and the election treatment is appointed
  baselines_list <- list()
  baselines_list$g.judge_sentencing <- "2"
  baselines_list$a.election_treatment <- "3"

#########################
# Support for the judge #
#########################
############## 
# Figure D.1 #
############## 
  # Variables to go in the model, ensuring they are alphabetically ordered in way we wish to present them
  survey2023$a.election_treatment <- survey2023$election_treatment
  survey2023$b.judge_name <- survey2023$judge_name
  survey2023$c.judge_race <- survey2023$judge_race
  survey2023$d.judge_margin <- survey2023$judge_margin
  survey2023$e.judge_tenure <- survey2023$judge_tenure
  survey2023$f.judge_party <- survey2023$judge_party
  survey2023$g.judge_sentencing <- survey2023$judge_sentencing

# Subset dataset to respondents who answered the support question 
  survey.support <- survey2023[!(is.na(survey2023$judge_support)),]
  survey.support$judge_support <- as.numeric(as.character(survey.support$judge_support))
  
# Estimate the basic AMCEs, all election types
  results.support.all.types <- amce(judge_support ~ a.election_treatment + b.judge_name  + c.judge_race + 
                                           d.judge_margin + e.judge_tenure + f.judge_party + g.judge_sentencing, 
                                         data=survey.support, baselines=baselines_list, cluster=F)

# Create figure

  plot(results.support.all.types, attribute.names = attribute.names.full,
                           label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1.25),
                           color="black", ylim=c(-1,1))
  
  ggsave("./figures_tables/cjoint-support-summer-2023.eps", width=7, height=5,device=cairo_ps)

#####################################################################
# In-text discussion of binary elections treatment from this survey #
#####################################################################
# Estimate the basic AMCEs, pooling election types together (HC2 matches amce output)
  results.support.binary <- lm(judge_support ~ judge_name + judge_race + I(election_treatment_binary==1) + judge_margin + judge_tenure
                               + judge_party + judge_sentencing, 
                               data=survey.support)
  results.support.binary <- coeftest(results.support.binary,vcovHC(results.support.binary,type="HC2"))
  # Estimate and p-value
  cat(round(results.support.binary["I(election_treatment_binary == 1)TRUE","Estimate"],2), sep = '\n', file = "./figures_tables/binary_coef_2023.tex")
  cat(p_ceiling(results.support.binary["I(election_treatment_binary == 1)TRUE","Pr(>|t|)"]), sep = '\n', file = "./figures_tables/binary_p_2023.tex")
    
##############################################
# Support for the judge as function of party #
##############################################
############## 
# Figure D.2 #
############## 
# Estimate the basic AMCEs, varying by respondent party
results.support.by.party <- amce(judge_support ~  a.election_treatment*Party + b.judge_name*Party  + c.judge_race*Party + 
                                   d.judge_margin*Party + e.judge_tenure*Party + f.judge_party*Party + g.judge_sentencing*Party, 
                                 data=survey.support, respondent.varying="Party", baselines=baselines_list, cluster=F,
                                 na.ignore = T)

  # Create figure
  figure.support.by.party <- plot(results.support.by.party, plot.display="interaction", attribute.names = attribute.names.full,
                                  label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1.25),
                                  color="black")
  
  ggsave("./figures_tables/cjoint-support-by-party-summer-2023.eps", width=7, height=5,device=cairo_ps)
  
##################################################
# Support for the judge as function of knowledge #
##################################################
############## 
# Figure D.3 #
##############
  # Estimate the basic AMCEs, varying by respondent knowledge, omitting respondents from NY due to ambiguity with "supreme court" question used to measure knowledge
  
  survey.support$`Knowledge Level` <- factor(survey.support$`Knowledge Level`,levels=c("Not Knowledgeable","Knowledgeable"))
  
  results.support.by.smart <- amce(judge_support ~ a.election_treatment*`Knowledge Level` + b.judge_name*`Knowledge Level` + c.judge_race*`Knowledge Level` + 
                                     d.judge_margin*`Knowledge Level` + e.judge_tenure*`Knowledge Level` + f.judge_party*`Knowledge Level` + g.judge_sentencing*`Knowledge Level`,na.ignore = T,
                                   data=survey.support[survey.support$STATE!="NY",], respondent.varying="Knowledge Level", baselines=baselines_list, cluster=F)

  # Create figure
  figure.support.by.smart <- plot(results.support.by.smart, plot.display="interaction", attribute.names = attribute.names.full,
                                  label.baseline=F, level.names=levels.test.full, main = "", xlab="Estimated effect", xlim=c(-1.25,1.25),
                                  color="black")

  ggsave("./figures_tables/cjoint-support-by-knowledge-summer-2023.eps", width=7, height=5,device=cairo_ps)
  
################################ 
# Tests of possible mechanisms #
################################
  # Binary variable for election treatment
  survey.support$got_election_treatment_binary <- ifelse(survey.support$election_treatment_binary == 1, 1, 0)
  # Measure of whether respondent got copartisan treatment, outpartisan treatment, or a no party listed treatment/respondent is an independent
  survey.support$copartisan <- NA
  survey.support[survey.support$judge_party==3 | survey.support$Party=="Independent" & !is.na(survey.support$Party),"copartisan"] <- "No Party Listed/Independent"
  survey.support[survey.support$judge_party==1 & survey.support$Party=="Republican" & !is.na(survey.support$Party),"copartisan"] <- "Copartisan"
  survey.support[survey.support$judge_party==2 & survey.support$Party=="Republican" & !is.na(survey.support$Party),"copartisan"] <- "Outpartisan"
  survey.support[survey.support$judge_party==1 & survey.support$Party=="Democrat" & !is.na(survey.support$Party),"copartisan"] <- "Outpartisan"
  survey.support[survey.support$judge_party==2 & survey.support$Party=="Democrat" & !is.na(survey.support$Party),"copartisan"] <- "Copartisan"
  survey.support$copartisan <- factor(survey.support$copartisan)

####################
####################
##### Efficacy #####
####################
####################
####################################
# External -- Do Politicians Care? #
####################################
#####################
# Figure 4, Panel a #
#####################
# Higher values indicate greater efficacy (so, they *disagree* that public officials don't care)
# Using interflex -- treating the efficacy response (moderator) as linear, interact with our binary treatment, determine marginal effect across range of moderator
  off_care_inter <- interflex(estimator="linear",
                               data=as.data.frame(survey.support),
                               Y="judge_support",
                               D="got_election_treatment_binary",
                               X="officials_dont_care",
                               Z=c("judge_name","judge_race","judge_margin","judge_tenure","judge_party","judge_sentencing"),
                               na.rm=T,full.moderate = F)
  # Plotting
  plot(off_care_inter,
       ylab="Effect of Election Treatment",
       xlab='"Public officials don\'t care what people like me think" \n (1=Strongly Agree, 5 = Strongly Disagree)')
  # Saving the plot
  ggsave("./figures_tables/officials_care_inter.eps",width=8,height=5,device=cairo_ps)

  # Getting estimate of treatment effect (to discuss in text)
  off_care_inter$Avg.estimate$`1`
  
##############
# Figure D.4 #
##############
# OLS, treat each level of efficacy as discrete category (factor) and interact with treatment
  do.officials.care.reg <- lm(judge_support ~ factor(officials_dont_care)*got_election_treatment_binary
                                                                            -got_election_treatment_binary+
                                                                            (judge_name + judge_race +
                                                                            judge_margin + judge_tenure 
                                                                          + judge_party + judge_sentencing), data=survey.support)
  
  # HC2 standard errors
  offic_care_fp <- coeftest(do.officials.care.reg,vcovHC(do.officials.care.reg,type="HC2"))
  # Efficacy * election treatment estimates
  offic_care_fp <- as.data.frame(offic_care_fp[grep("election_treatment",rownames(offic_care_fp)),])
  # Levels of efficacy measure
  offic_care_fp$moderator <- gsub("[^0-9]","",rownames(offic_care_fp))
  # Plotting estimated effects
  ggplot(offic_care_fp,aes(x=moderator,y=Estimate,
                           ymax=Estimate+1.96*`Std. Error`,
                           ymin=Estimate-1.96*`Std. Error`))+
    geom_hline(yintercept=0,colour="grey50",size=1.25,linetype=2)+
    geom_point(size=3)+
    geom_linerange(size=1.25)+
    ylab("Effect of Election Treatment")+
    xlab('"Public officials don\'t care what people like me think" \n (1=Strongly Agree, 5 = Strongly Disagree)')+
    theme_minimal()
  # Saving the plot
  ggsave("./figures_tables/officials_care_discrete.eps",width=8,height=5,device=cairo_ps)

#######################################
# Internal -- Politics is Complicated #
#######################################
#####################
# Figure 4, Panel b #
#####################
# Coded so higher values indicate greater efficacy (so, they *disagree* that they can't understand politics)
  understand_inter <- interflex(estimator="linear",
            data=as.data.frame(survey.support),
            Y="judge_support",
            D="got_election_treatment_binary",
            X="cant_understand_politics",
            Z=c("judge_name","judge_race","judge_margin","judge_party","judge_tenure","judge_sentencing"),
            na.rm=T,full.moderate = F)
  # Getting p-value on estimate of slope (to discuss in text)
  cat(as.character(format(nsmall=2,p_ceiling(summary(understand_inter$model.linear)$coefficients["DX.Group.2","Pr(>|t|)"]))),
      sep = '\n', file = "./figures_tables/internal_p.tex")
  # Plotting
  plot(understand_inter,
       ylab="Effect of Election Treatment",
       xlab='"Sometimes politics and government seem so complicated that \n a person like me can\'t really understand what\'s going on" \n (1=Strongly Agree, 5 = Strongly Disagree)')
  # Saving the plot
  ggsave("./figures_tables/understand_inter.eps",width=8,height=5,device=cairo_ps)

##############
# Figure D.5 #
##############
# OLS
  understand_reg <- lm(judge_support ~ factor(cant_understand_politics)*got_election_treatment_binary
                              -got_election_treatment_binary+
                                (judge_name + judge_race +
                                   judge_margin + judge_tenure 
                                 + judge_party + judge_sentencing), data=survey.support)
  # HC2 standard errors
  understand_fp <- coeftest(understand_reg,vcovHC(understand_reg,type="HC2"))
  # Efficacy * election treatment estimates
  understand_fp <- as.data.frame(understand_fp[grep("election_treatment",rownames(understand_fp)),])
  # Levels of efficacy measure
  understand_fp$moderator <- gsub("[^0-9]","",rownames(understand_fp))
  # Plotting estimated effects
  ggplot(understand_fp,aes(x=moderator,y=Estimate,
                           ymax=Estimate+1.96*`Std. Error`,
                           ymin=Estimate-1.96*`Std. Error`))+
    geom_hline(yintercept=0,colour="grey50",size=1.25,linetype=2)+
    geom_point(size=3)+
    geom_linerange(size=1.25)+
    ylab("Effect of Election Treatment")+
    xlab('"Sometimes politics and government seem so complicated that \n a person like me can\'t really understand what\'s going on" \n (1=Strongly Agree, 5 = Strongly Disagree)')+
    theme_minimal()
  # Saving the plot
  ggsave("./figures_tables/understand_discrete.eps",width=8,height=5,device=cairo_ps)

###################################
###################################
###### Support for democracy ######
###################################
###################################
################################
# Local Elections are Worth It #
################################
#####################
# Figure 5, Panel b #
#####################
# Coded so higher values indicate greater democratic support (so, they *disagree* that local elections aren't worth bothering with)
  worth_inter <- interflex(estimator="linear",
                                data=as.data.frame(survey.support),
                                Y="judge_support",
                                D="got_election_treatment_binary",
                                X="local_elections_not_worth",
                                Z=c("judge_name","judge_race","judge_margin","judge_tenure","judge_party","judge_sentencing"),
                                na.rm=T,full.moderate = F)
  # Plotting
  plot(worth_inter,
       ylab="Effect of Election Treatment",
       xlab='"Many local elections aren\'t important enough to bother with." \n (1=Strongly Agree, 5 = Strongly Disagree)')
  # Saving the plot
  ggsave("./figures_tables/worth_inter.eps",width=8,height=5,device=cairo_ps)

##############
# Figure D.7 #
##############  
# OLS
  worth_reg <- lm(judge_support ~ factor(local_elections_not_worth)*got_election_treatment_binary
                       -got_election_treatment_binary+
                         (judge_name + judge_race +
                            judge_margin + judge_tenure 
                          + judge_party + judge_sentencing), data=survey.support)
  # HC2 standard errors, democracy * election treatment estimates, levels of democracy measure
  worth_fp <- coeftest(worth_reg,vcovHC(worth_reg,type="HC2"))
  worth_fp <- as.data.frame(worth_fp[grep("election_treatment",rownames(worth_fp)),])
  worth_fp$moderator <- gsub("[^0-9]","",rownames(worth_fp))
  # Plotting estimated effects
  ggplot(worth_fp,aes(x=moderator,y=Estimate,
                           ymax=Estimate+1.96*`Std. Error`,
                           ymin=Estimate-1.96*`Std. Error`))+
    geom_hline(yintercept=0,colour="grey50",size=1.25,linetype=2)+
    geom_point(size=3)+
    geom_linerange(size=1.25)+
    ylab("Effect of Election Treatment")+
    xlab('"Many local elections aren\'t important enough to bother with." \n (1=Strongly Agree, 5 = Strongly Disagree)')+
    theme_minimal()
  # Saving the plot
  ggsave("./figures_tables/worth_discrete.eps",width=8,height=5,device=cairo_ps)

######################################
# Elections Force Government to Care #
######################################
#####################
# Figure 5, Panel a #
#####################
# Coded so higher values indicate greater democratic support (so, they *agree* that elections make the government care)
  govtcare_inter <- interflex(estimator="linear",
                           data=as.data.frame(survey.support),
                           Y="judge_support",
                           D="got_election_treatment_binary",
                           X="elections_make_govt_care",
                           Z=c("judge_name","judge_race","judge_margin","judge_tenure","judge_party","judge_sentencing"),
                           na.rm=T,full.moderate = F)
  # Plotting
  plot(govtcare_inter,
       ylab="Effect of Election Treatment",
       xlab='"Elections make the government pay attention to what the people think." \n (1=Strongly Disagree, 5 = Strongly Agree)')
  # Saving the plot
  ggsave("./figures_tables/care_inter.eps",width=8,height=5,device=cairo_ps)
 
##############
# Figure D.6 #
############## 
# OLS
  govtcare_reg <- lm(judge_support ~ factor(elections_make_govt_care)*got_election_treatment_binary
                  -got_election_treatment_binary+
                    (judge_name + judge_race +
                       judge_margin + judge_tenure 
                     + judge_party + judge_sentencing), data=survey.support)
  # HC2 standard errors, democracy * election treatment estimates, levels of democracy measure
  govtcare_fp <- coeftest(govtcare_reg,vcovHC(govtcare_reg,type="HC2"))
  govtcare_fp <- as.data.frame(govtcare_fp[grep("election_treatment",rownames(govtcare_fp)),])
  govtcare_fp$moderator <- gsub("[^0-9]","",rownames(govtcare_fp))
  # Plotting estimated effects
  ggplot(govtcare_fp,aes(x=moderator,y=Estimate,
                      ymax=Estimate+1.96*`Std. Error`,
                      ymin=Estimate-1.96*`Std. Error`))+
    geom_hline(yintercept=0,colour="grey50",size=1.25,linetype=2)+
    geom_point(size=3)+
    geom_linerange(size=1.25)+
    ylab("Effect of Election Treatment")+
    xlab('"Elections make the government pay attention to what the people think." \n (1=Strongly Disagree, 5 = Strongly Agree)')+
    theme_minimal()
  # Saving the plot
  ggsave("./figures_tables/care_discrete.eps",width=8,height=5,device=cairo_ps)
  
###############################
# Satisfaction with Democracy #
###############################
#####################
# Figure 5, Panel c #
#####################
# Coded so higher values indicate greater democratic support (so, they *are satisfied* with democracy)
  democracy_inter <- interflex(estimator="linear",
                            data=as.data.frame(survey.support),
                            Y="judge_support",
                            D="got_election_treatment_binary",
                            X="democracy_satisfaction",
                            Z=c("judge_name","judge_race","judge_margin","judge_tenure","judge_party","judge_sentencing"),
                            na.rm=T,full.moderate = F)
  # Getting p-value on estimate of slope (to discuss in text)
  cat(as.character(format(nsmall=2,p_ceiling(summary(democracy_inter$model.linear)$coefficients["DX.Group.2","Pr(>|t|)"]))),
      sep = '\n', file = "./figures_tables/satisfaction_p.tex")
  # Plotting
  plot(democracy_inter,
     ylab="Effect of Election Treatment",
     xlab='"How satisfied are you with the way democracy is working in the United States?" \n (1=Very Dissatisfied, 5 = Very Satisfied)')
  # Saving the plot
  ggsave("./figures_tables/democracy_inter.eps",width=8,height=5,device=cairo_ps)

##############
# Figure D.8 #
##############
# OLS
  satisfaction_reg <- lm(judge_support ~ factor(democracy_satisfaction)*got_election_treatment_binary
                         -got_election_treatment_binary+
                           (judge_name + judge_race +
                              judge_margin + judge_tenure
                            + judge_party + judge_sentencing), data=survey.support)
  # HC2 standard errors, democracy * election treatment estimates, levels of democracy measure
  satisfaction_fp <- coeftest(satisfaction_reg,vcovHC(satisfaction_reg,type="HC2"))
  satisfaction_fp <- as.data.frame(satisfaction_fp[grep("election_treatment",rownames(satisfaction_fp)),])
  satisfaction_fp$moderator <- gsub("[^0-9]","",rownames(satisfaction_fp))
  # Plotting estimated effects
  ggplot(satisfaction_fp,aes(x=moderator,y=Estimate,
                       ymax=Estimate+1.96*`Std. Error`,
                       ymin=Estimate-1.96*`Std. Error`))+
  geom_hline(yintercept=0,colour="grey50",size=1.25,linetype=2)+
  geom_point(size=3)+
  geom_linerange(size=1.25)+
  ylab("Effect of Election Treatment")+
  xlab('"How satisfied are you with the way democracy is working in the United States?" \n (1=Very Dissatisfied, 5 = Very Satisfied)')+
  theme_minimal()
  # Saving the plot
  ggsave("./figures_tables/democracy_discrete.eps",width=8,height=5,device=cairo_ps)

#################################################################################
# Putting the efficacy and support for democracy results into a stargazer table #
#################################################################################
#############
# Table D.2 #
#############
# Efficacy (external)
survey.support$moderator <- survey.support$officials_dont_care
officials_care_lm <- lm(judge_support ~ moderator*got_election_treatment_binary +
                                      judge_name + judge_race + judge_margin + judge_tenure + judge_party + judge_sentencing, 
                                      data=survey.support)
# Efficacy (internal)
survey.support$moderator <- survey.support$cant_understand_politics
understand_lm <- lm(judge_support ~ moderator*got_election_treatment_binary +
                                      judge_name + judge_race + judge_margin + judge_tenure + judge_party + judge_sentencing, 
                                      data=survey.support)
# Democracy (local elections)
survey.support$moderator <- survey.support$local_elections_not_worth
worth_lm <- lm(judge_support ~ moderator*got_election_treatment_binary +
                                      judge_name + judge_race + judge_margin + judge_tenure + judge_party + judge_sentencing, 
                                      data=survey.support)
# Democracy (elections force)
survey.support$moderator <- survey.support$elections_make_govt_care
govtcare_lm <- lm(judge_support ~ moderator*got_election_treatment_binary +
                                     judge_name + judge_race + judge_margin + judge_tenure + judge_party + judge_sentencing, 
                                      data=survey.support)
# Democracy (satisfaction)
survey.support$moderator <- survey.support$democracy_satisfaction
democracy_lm <- lm(judge_support ~ moderator*got_election_treatment_binary +
                                     judge_name + judge_race + judge_margin + judge_tenure + judge_party + judge_sentencing, 
                                      data=survey.support)
# List of all regressions
moderator_regs_list <- list(officials_care_lm,
                            understand_lm,
                            worth_lm,
                            govtcare_lm,
                            democracy_lm)
# Stargazer table
figure_mods <- stargazer(moderator_regs_list,
                         se=lapply(list(officials_care_lm,understand_lm,worth_lm,govtcare_lm,democracy_lm),function(x) sqrt(diag(vcovHC(x,type="HC2")))),
                         dep.var.labels = "Support for Judge",
                         column.labels = c("External Efficacy","Internal Efficacy","Local Elections","Elections Force","Democracy Support"),
                             covariate.labels = c("Moderator",
                                                  "Elected",
                                                  "Woman",
                                                  "White",
                                                  "Hispanic",
                                                  "60\\% Margin",
                                                  "80\\% Margin",
                                                  "5 Year Tenure",
                                                  "15 Year Tenure",
                                                  "Democrat",
                                                  "No Party Listed",
                                                  "6 Year Avg. Sentence",
                                                  "9 Year Avg. Sentence",
                                                  "Moderator $\\times$ Elected",
                                                  "Constant"),
                             keep.stat=c("n","adj.rsq"),
                             notes.append = FALSE,
                             table.layout ="-ldc-t-as-n",font.size="footnotesize",
                             star.char="*",star.cutoffs=c(0.1,0.05),notes.label = "",
                             column.sep.width="-0.20cm",no.space=T,
                             notes="\\parbox[t]{1\\textwidth}{\\footnotesize \\textit{Note}: Estimates are from OLS regressions with robust standard errors in parentheses. 
                             The unit of observation is the respondent. $^{*}p<0.1$, $^{*}p<0.05$.}",
                             title="Political Efficacy and Attitudes Toward Democracy as Moderators",
                             label="efficacy_regs"
)
# Saving the table
cat(figure_mods, sep = '\n', file = "./figures_tables/efficacy_table.tex")

#############################################
#############################################
###### Perceived ideological proximity ######
#############################################
#############################################
###########
# Table 4 #
###########
# Punitiveness 
# OLS, control for punitiveness distance
punitiveness.distance.reg <- lm(judge_support ~ judge_name + judge_race 
                              + got_election_treatment_binary + 
                                judge_margin + judge_tenure 
                              + copartisan + judge_sentencing
                              + punitiveness_distance, data=survey.support)

# Ideology 
# OLS, ideological distance as outcome
perc_ideology.distance.5pt.reg <- lm(ideology_distance_5 ~ judge_name + judge_race 
                                + got_election_treatment_binary + 
                                  judge_margin + judge_tenure 
                                + copartisan + judge_sentencing, data=survey.support)

# OLS, control for ideological distance
ideology.distance.5pt.reg <- lm(judge_support ~ judge_name + judge_race 
                         + got_election_treatment_binary + 
                           judge_margin + judge_tenure 
                         + copartisan + judge_sentencing
                         + ideology_distance_5, data=survey.support)

# Putting the proximity results into stargazer tables
# Lists of outcomes
proximity_perc_regs <- list(perc_ideology.distance.5pt.reg)
proximity.regs <- list(punitiveness.distance.reg, ideology.distance.5pt.reg)
# Stargazer table
proximity_output <- stargazer(proximity_perc_regs,proximity.regs,
                              se=lapply(list(perc_ideology.distance.5pt.reg,
                                             punitiveness.distance.reg,
                                             ideology.distance.5pt.reg),function(x) sqrt(diag(vcovHC(x,type="HC2")))),
                              dep.var.labels = c(#"Punitive Proximity",
                                                 "Ideological Distance","Judge Support"),
                              order=c(4,1,2,3,5:14),
                             covariate.labels = c("Elected",
                                                  "Woman",
                                                  "White",
                                                  "Hispanic",
                                                  "60\\% Margin",
                                                  "80\\% Margin",
                                                  "5 Year Tenure",
                                                  "15 Year Tenure",
                                                  "No Party Listed/Independent",
                                                  "Outpartisan",
                                                  "6 Year Avg. Sentence",
                                                  "9 Year Avg. Sentence",
                                                  "Punitive Distance (0-2)",
                                                  "Ideological Distance (0-4)",
                                                  "Constant"),
                              keep.stat=c("n","adj.rsq"),
                             notes.append = FALSE,
                             table.layout ="-ld-t-as-n",font.size="footnotesize",
                             star.char="*",star.cutoffs=c(0.1,0.05),notes.label = "",
                             column.sep.width="0cm",no.space=T,
                             notes="\\parbox[t]{0.7\\textwidth}{\\footnotesize \\textit{Note}: Estimates are from OLS regressions with robust standard errors in parentheses.
                             The unit of observation is the respondent.
                             Punitive distance is a 3-point measure (0-2); ideological distance is a 5-point measure (0-4). 
                             Higher values indicate greater perceived distance. $^{*}p<0.1$, $^{*}p<0.05$.}",
                             title="Perceived Political Proximity and The Value of Judicial Elections",
                             label="proximity_regs"
)
# Saving the table
cat(proximity_output, sep = '\n', file = "./figures_tables/proximity_table.tex")

#####################################
# Elections and shared partisanship #
#####################################
#######################
# Figure E.1, Panel b #
#######################
# Subset to respondents who are partisans and got a partisan judge treatment
survey.support.partisans <- survey.support[which(survey.support$Party %in% c("Democrat","Republican") & survey.support$judge_party!=3),]
# Copartisan variable
survey.support.partisans$judge_copartisan <- ifelse((survey.support.partisans$judge_party == 2 & survey.support.partisans$Party == "Democrat") |
                                                      (survey.support.partisans$judge_party == 1 & survey.support.partisans$Party == "Republican"), 1, 0)

# OLS, support on indicator for the binary election treatment interacted with copartisanship
f2.a.reg <- lm(judge_support ~ judge_copartisan*election_treatment_binary, data=survey.support.partisans)
# Predicted values
f2.a.dataset <- data.frame(election_treatment_binary=as.factor(c(1,2,1,2)), judge_copartisan=c(0,0,1,1))
f2.a.fitted.robust <- predict.rob(x=f2.a.reg,newdata=f2.a.dataset,clcov=vcovHC(f2.a.reg,type="HC2"))

# Creating the plot
coef <- c(f2.a.fitted.robust$fit[1],f2.a.fitted.robust$fit[3],f2.a.fitted.robust$fit[2],f2.a.fitted.robust$fit[4])
se <- c(f2.a.fitted.robust$se.fit[1],f2.a.fitted.robust$se.fit[3],f2.a.fitted.robust$se.fit[2],f2.a.fitted.robust$se.fit[4])
lb <- coef-1.96*se
ub <- coef+1.96*se
var <- c("Elected","Elected","Appointed","Appointed")
party <- c("Outpartisan","Copartisan","Outpartisan","Copartisan")

results <- data.frame(coef,se,lb,ub,var,party)
results$var <- reorder(results$var,rep(1:2,each=2))

ggplot(results, aes(x=coef,y=var,color=party,group=party)) +
  geom_point(size=4) + 
  geom_segment(aes(x=lb,xend=ub,y=var,yend=var),size=1.25)+
  scale_color_grey(end=0.6)+
  theme_bw()+
  xlab("Predicted Support")+
  labs(colour="")+
  theme(legend.position="bottom",
        axis.title.y = element_blank())
# Saving the plot
ggsave("./figures_tables/elections-moderate-differences-study-two.eps",width=6,height=4,device=cairo_ps)
